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ABSTRACT 


In  fitting  a  transfer  function  model  to  input/output  data,  it  is  sometimes 
possible  to  select  a  particular  input  with  the  purpose  of  minimizing  estimation 
error.  In  this  report,  some  aspects  of  this  optimal  input  synthesis  problem 
are  examined  for  a  class  of  linear  multistage  systems.  In  particular,  an 
input  is  selected  from  a  class  of  equal- energy  inputs  which  maximizes  the 
trace  of  the  Fisher  information  matrix  with  respect  to  the  input  sensitivity 
parameters,  where  it  is  assumed  that  the  transient  response  and  prediction 
accuracy  of  the  unforced  system  is  known.  The  optimal  input  is  compared  to 
a  purely  random  input  which  has  the  same  average  energy,  and  numerical 
results  for  a  second-order  system  are  provided,  input  synthesis  with  regard 
to  the  other  model  parameters  which  determine  its  unforced  behavior  is 
hampered  by  the  fact  that  the  Fisher  information  matrix  is  a  function  of  these 
parameters. 
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NOMENCLATURE 


The  following  is  a  list  of  symbols  used  in  the  report.  Underlined  letters 
denote  vectors. 


h 

I 

n 

J 

k 

L{e) 

L'(e) 

M 

N 

n 

N  (m,Q) 
n  — 


Description 


Matrix  of  input  coefficients 
Partial  of  A  with  respect  to 
Matrix  of  noise  input  coefficients 

Covariance  of  the  steady -state  measurement  residuals 
Sample  covariance  of  the  unconstrained  likelihood  function 
N  by  N  symmetric  matrix  which  is  a  sum  of  similarity 
transforms  of  S 

Measurement  matrix  (a  vector  for  scalar  measurements) 
Identity  matrix  of  dimension  n  by  n 
Cost  function 

Steady- state  Kalman  filter  gain 
Unconstrained  likelihood  function 
Constrained  likelihood  function 
The  Fisher  information  matrix 
Total  number  of  measurements 
Dimension  of  x(i) 

Muxtivariant  normal  distribution  of  dimension  n  with 
mean  m  and  covariance  Q 


VI 
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Symbol 

Q 

r(i) 

u 

u(i) 

u-{i) 

S(i) 

V 

v(i) 

w(i) 

x(i) 

x'(i) 


z(i) 

a 

1 

c 

e 

^ave 


NOMENCLATURE  (Continued) 

t 

>  i 

Description 

Covariance  of  the  random  input  r{i) 

! 

Random  input  at  index  i 
Composite  vector  of  inputs  ' 

Input  at  index  i  > 

Optimal  input  at  index  i 

Random  input  at  index  i  i 

Composite  vector  of  measurement' residuals 
Measurement  residual  at  index  i 
Observation  vector  at  index  i 
State  vector  at  index  i 

•  I 

One -step -ahead  predicted  estimate  of  the  state  vector 
at  index  i  i 

I 

Composite  vector  of  measurements 

Output  at  index  i  , 

Total  energy  of  input 
Input  sensitivity  vector 

Efficiency  of  the  random  input  compared  to  optimal  input 

Unknown  parameter  vector 

Lagrange  multiplier 

Average  of  eigenvalues  of  F  ' 


I 
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NOMENCLATURE  (Continued) 

1  Description 

Maxiinum  eigenvalue  of  F 

Covariance 'of  unconstrained  likelihood  function 

State  transition  matrix 

Parameters  in  state  transition  matrix 


I.  INTRODUCTION 


Measurements  of  the  input  and  output  of  an  isolated  physical  system 
are  made  at  uniformly  spaced  instants  of  time.  The  sets  of  N  measurements 
of  the  input  and  output  are  denoted  {u(i):  i  -  1,  •  •  •  ,  N)  and  {z(i);  1=1, 

•  •  •  ,  N},  respectively.  The  measurements  are  assumed  to  be  error-free. 
Furthermore,  {z(i):  i  =  1,  •  •  •  ,  N)  are  assumed  to  be  the  output  of  a 
multistage  time-invariant  linear  system  forced  by  the  inputs  {u(i):  i  =  1, 

•  •  •  ,  N}  ana  a  postulated  sequence  of  purely  random  inputs  which  are  assumed 
to  arise  due  to  modeling  errors.  The  multistage  model  contains  unknown 
parameters  which  must  be  identified.  The  general  problem  of  estimating  the 
unknown  parameters  of  this  linear  stochastic  model  from  the  input/output 
records  is  described  in  ( l] . 

In  some  applications  it  is  possible  to  apply  a  particular  input  sequence 
{u(i):  i  =  1,  •  •  ,  N}  specifically  for  the  purpose  of  enhancing  the  estimates 

of  the  unknown  model  parameters^.  The  basic  conditions  are: 

1)  The  system  to  be  identified  is  out  of  service:  that  is, 
one  has  complete  control  over  the  input  during  the 
identification  phase. 

2)  The  input  resource  (energy)  is  scarce. 

If  there  is  no  penalty  or  constraint  on  the  input,  then  optima’  input 
synthesis  is  unnecessary.  This  is  because  any  input  which  excites  all  the 

^This  problem  is  the  dynamic  analog  to  experimental  design  in  statistics  [Z], 
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modes  of  the  system,  if  applied  long  enough,  is  generally  sufficient  to  allow 
convergence  of  the  identification.  There  are  several  systems  for  which 
there  is  a  time,  magnitude,  or  energy  constraint  on  the  input.  One  example 
is  a  structure  where  the  input  is  a  thrust  profile  from  a  rocket  engine  which 
has  a  finite  supply  of  fuel.  Another  example  is  process  control  where  the 
amount  of  time  spent  in  the  parameter  identification  phase  can  be  very  costly. 

The  problem  of  selecting  input  signals  to  enhance  the  system  identifi¬ 
cation  process  has  been  examined  by  several  authors.  Aoki  and  Staley  [3] 
have  obtained  some  general  results  for  a  wide  class  of  discrete-time  systems 
as  well  as  justifying  the  use  of  the  trace  of  the  Fisher  information  matrix  as 
a  performance  index.  This  paper  retains  the  use  of  this  convenient  per¬ 
formance  index,  but  formulates  the  problem  differently  by  using  the  Kalman 
filter  representation  of  th  system.  Levadi  [4]  considers  identifying  non¬ 
stationary,  continuous  systems  with  optimum  energy- constrained  and  time- 
constrained  input  signals.  The  only  source  of  noise  considered,  however,  is 
output  observation  noise.  Gagliardi  [s]  attacks  the  problem  in  a  completely 
different  manner  by  using  multiple  hypothesis  testing. 

In  this  paper  a  Kalman  filter  representation  results  in  a  considerable 
simplification  of  the  problem.  It  is  assumed  that  the  input  sequence  is 
selected  from  a  class  of  possible  input  sequences  with  equal  total  energy  or 
equal  total  average  energy  if  the  input  is  random.  This  assumption  attempts 
to  bound  the  output  of  the  system  to  prevent  forcing  it  into  a  region  where  it 
is  no  longer  linear.  The  optimal  input  is  the  sequence  in  this  class  which 
minimizes  the  trace  of  the  Fisher  information  matrix.  Some  investigators  ha\ 


suggested  the  use  of  a  purely  random  input  to  drive  the  system  during  the 
identification  phase  [6]  .  A  comparative  efficiency  between  a  purely  random 
input  and  the  optimal  input  for  the  class  of  inputs  studied  is  defined,  and  some 
numerical  results  for  a  second-order  model  are  presented. 

It  should  be  emphasized  at  the  outset  that  there  is  a  fundamental  limita¬ 
tion  on  the  solution  of  the  general  optimal  input  synthesis  problem:  the  Fisher 
inforr.'.  cion  matrix  is  generally  a  function  of  the  unknown  parameters.  Hence, 
one  is  faced  with  the  impossible  situation  of  requiring  the  true  values  of  the 
unknown  parameters  in  order  to  find  the  optimal  input.  In  the  general  case, 
therefore,  any  physically  realizable  efficient  input  synthesis  technique  must 
be  adaptive  --as  better  estimates  of  the  parameters  are  obtained,  the  input 
is  made  to  approach  the  appropriate  optimal  input.  There  is  an  exception  in 
linear  systems  which  is  examined  in  this  report:  the  Fisher  inforniation 
matrix  is  not  a  function  of  the  input  parameters. 
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n.  THE  MODEL 


A  linear  multistage  model  which  relates  the  input  {u(i): 
to  the  output  {z(i):  i  =  1,  .  . . ,  N}  is 

x(i  +  1)  =  $x(i)  +  Y  u(i)  +  r^(i) 
z(i)  =  h^  x(i) 


i  =  1 ,  .  .  .  ,  N) 


(1) 

(2) 


where  {^(i):  i  =  1,  . .  . ,  N)  is  a  set  of  independent  vector -valued  random 
'•ariables  each  of  which  are  distributed  N^(0^,  Q),  whore  Q  is  a  diagonal 
matrix  with  non-negative  elements.  The  initial  value  x(i)  is  assumed  to 
be  0.  A  canonical  form  of  Eqs.  (1)  and  (2)  utilizes  the  matrix  forms 
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ro  1 1  ,1 

_ 

-  ,  n-l 

h  = 

1 

1 

1 

1 

.  0  . 

The  unmeasured  states  {xj(i):  i  =  1 ,  .  ,  ,  ,  N;  j  =  2,  .  .  .  ,  n)  can  be  estimated  by 
the  corresponding  steady-state  Kalman  one -step-ahead  predictor  [7]  given 
by 


x'(i  +  1)  =  $x'(i)  $kv(i)  +  Y  u(i ) 

(3) 

z(i )  =  h"^  x'(i)  +  v(i) 

(•!) 
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1 

k  is  the  steady-state  Kalmain  filter  gain.  It  is  straightforward  to  estimate 

the  unknown  model  parameter^  and  k  from  the  'filter  residuals  [l].  If  one  , 

assumes  steady-state  , conditions,  the  measurement  residuals  {v(i): 

i  =  1, . . . ,  N'}  are  iUidependent  random  variables  each  of  which  is  distributed 

N(0,  b)  as  a  result  of  the  innovations iproperty  of  the  Kalmkn  filter  [8], 

.  1  *  ** 
Equations  (3)  and  (4)  can  also  be  u'sed  to  represent  the  relationship  between 

the  input  {u(i/:  i  =  1, . . . ,  N}'and  the  output  {z(i)t  i  =  1,’. . .  ^  N},  since  they  i 

are  statistically  equivalent  [7].  ‘ 


I  ! 

I 

I  ' 
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m.  DERIVATION  OF  AN  UNCONSTRAINED  LIKELIHOOD  FUNCTION 


•  ,  I  I 

The  likelihood  function' which  is  most  convenient  for  computing  the 
maximum -likelihood  estimates  of  the  unknown  model  parameters  is 


N 


■L'ie)  =  -  f  log.  b  v^i)' 


(5) 


i=l 


where  {v(i):  i  =  1*, .  .  .  ,  N}  are  constrained  to  satisfy  Eq.  (3),  and  where  G  is 
the  vector  of  the  unknown  model  parameters  [l].  The  constraint  on  the 

I 

measurement  residuals  is  cumbersome  to  work  with  analytically,  so  an 

equivalent  likelihood  function , is  derived  in  this  section  from  the  distribution 

,  » 

of  the  rAeasurements  {zH):  i  =  1, ...»  k},  wnich  do  not  have  auxiliary 
constraints.  ' 

First,  Eq,  (3)  is  solved  in  closed  form: 


(  n-1  n-1 ' 

■J  =  l  j  =  l. 


(6) 


Since  the  measurements  are  assumed  to  be  error-free,  it  follows  that 


z(n)  =_h^^  (x'(n)  +  k  \  {n)) 


(7) 
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Substitute  Eq.  (6)  into  Eq.  (7)  to  obtain  a  closed  form  expression  for  z(n); 

n-1  n 

z(n)  =  Y  u(j)  +  ^  ^  liv(j) 

j=l  j=l 

Define  the  following  composite  N  vectors: 


u#  :  : 


It  follows  from  Eq.  (8)  that 


z  =  A  u  +  B  V 


where  A  and  B  are  lower  triangular  N  by  N  matrices  defined  as 


;  f  £  m 


,  T  .f-1  -m  ^ 

h  $  y  •,  i  >  ni 


f  <  m 


:  f  =  m 


,  T  .i-m  ,  .  ^ 

h  4>  k  ;  f  >  m 
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Equation  (9)  shows  that  the  measurements  ^are  linear  combinations  of  the 
inputs  u  and  the  measurement  residuals  v  (which  are  assumed  to  result 
from  modeling  errors),  assuming  zero  initial  conditions. 

To  simplify  the  notation,  define  the  n  vector. 


,  .  A  /»T.n  , 
w(n)  =  (®  )  h 

and  substitute  it  into  Eqs.  (10)  and  (11): 

i  <  m 


^  >  m 


"im 


Ly  w  (i  -1  -m) 


(12) 


(13) 


w(i 


■m) 


f  <  m 
i  =  m 

i  >  m 


(H) 


From  the  innovations  property  v  is  distributed  Nj^(0^,  bl)  and  so  by  Eq.  (9) 

T 

^  is  distributed  Nj^(Aii,  bB  B  ).  The  likelihood  function  E(0^)  for  the 
multivariate  normal  distribution  is  a  standard  result  in  statistics  (2|: 

L(0)  =  -  j  logjSi  -  ™  Tr  [S"^C] 

A  T 

where  S  “  b  B  B 

C  ^  (^  -  A  u)  (^  -  A  u)^  (15) 
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IV.  AN  INPUT  SYNTHESIS  PROBLEM 


In  the  most  general  modeling  problems  the  values  of  k,  and  b 
are  uncertain.  For  the  remainder  of  this  section,  however,  it  is  assumed 
that  the  transient  response  of  the  system  (^)  and  the  one -step -ahead  pre¬ 
diction  accuracy  of  the  unforced  model  (^)  are  known  or  have  been  previously 
identified,  but  that  the  input  sensitivity  (^)  is  unknown.  The  true  value 
of  b  is  not  required  in  this  case.  It  is  also  assum  J  that  the  input  ti  is  to  be 
selected  from  a  class  of  inputs  such  that 

u^  11  =  O'  (16) 

where  ot  is  the  total  input  energy.  The  question  is:  what  form  of  the  con¬ 
strained  input  u  will  yield  the  most  "information  about  V'?  The  measure 
of  "information  about  Y_"  is  defined  as  the  trace  of  the  Fisher  information 
matrix  [3]. 

The  Fisher  information  matrix  is  a  measure  of  the  dispersion  of  the 
gradient  of  the  likelihood  function  v/ith  respect  to  the  unknown  parameters. 

At  one  extreme,  if  there  is  a  one-to-one  correspondence  between  the  values 
of  the  unknown  parameters  and  the  measurements,  then  the  information  is 
a  maximum.  At  the  other  extreme,  if  there  is  no  correspondence,  the 
information  matrix  is  zero.  See  [3]  for  additional  discussion. 

The  Fisher  information  mat  .'ix  (denoted  by  M)  for  a  multivariate  normal 
distribution  is  a  standard  result  in  statistics  [9]. 

Preceding  pege  blank 
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The  elements  of  M  are 
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To  proceed,  define  a  scalar  cost  fxinctional  J  which  is  proportional  to  Tr[Mj 
and  adjoin  Eq.  (17)  with  a  Lagrange  multiplier: 


J  =  F  u  +  \  (u^u  -  a) 


(21) 


The  gradient  of  J  is 


V  J  =  F  u  -  \  u  (22) 

_u  —  — 

Setting  the  gradient  to  zero,  the  result  is 

Fu  =  \u  (23) 

It  follows  from  Eq.  (23)  that  the  input  sequence  which  maximizes  J  is  the 
eigenvector  of  F  (denoted  ^  )  which  corresponds  to  the  largest  eigenvalue  of 
F  (denoted  \  ).  It  also  follows  from  Eq.  (23)  that 


max  Tr[M) 


max 


a 


(24) 


Note  that  neither  J  nor  u'  is  a  function  of  b. 


V.  COMPARISON  OF  THE  OPTIMAL  INPUT  WITH  AN 
EQUIVALENT  RANDOM  INPUT 


Some  investigators  have  used  random  inputs  to  drive  the  system  because 
a  random  input  is  most  likely  to  excite  all  the  modes  of  the  system  and  is 
very  easy  to  generate  [6],  In  this  section,  a  purely  random  input  is  compared 
(in  a  method  defined  below)  with  the  optimal  input  derived  in  the  last  section. 
First,  let  {u(i):  i  =  1,  .  .  .  ,  N}  denote  an  input  sequence  which  is  a  set  of 
independent  random  numbers  with  mean  zero,  such  that 

Efu^(i)]=§  (25) 


Therefore  u  has  an  average  energy  equal  to  a  ,  Furthermore,  suppose  that 
{u(i):  i  =  i,  .  .  .  ,  N)  are  independent  of  the  measurement  residuals 
(v(i):  i  =  1,  .  .  .  ,  N}.  The  efficiency  of  the  random  input  is  defined  as 
the  ratio 


I^m] 

max  Tr[M] 
u 


(26) 


where  the  overbar  denotes  the  expected  value  of  Tr[Mj  using  the  random  input. 
Since  u  is  assumed  to  be  independent  of  v,  it  follows  that 

E[Tr|Mj]  =  Tr[E[Ml]--iE[u’'Fu|  --  Tr[F)  = (271 


Preceding  page  blank 
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where  denotes  the  average  of  the  eigenvalues  of  F.  It  follows  from 

Eq.  (24)  that  the  efficiency  of  the  random  input  compared  ^to  the  optimal 

I 

\ 


input  is 


VI.  EXAMPLE 


I  I  ... 

■  The  purpose  of  this  example  is  two-fold:  first,  to  illustrate  the  theory- 
developed  in  this  report  and,  second,  to  try  tp  obtain  some  insight  into  the 

{  I 

form  of  the  loptimal  input  for  a  certajin  class  of  identification  problems.  The 
following  Kalman  filter  representation  of  a  second-order  system  is  examined: 


} 

0 

'  1 

1 

0 

1 

1 

x'(i  +  1) 

-^1. 

-h 

5'(i)  ■+ 

^2 

v(i)  + 

^2 

'  1_ 

I 


(29) 


To  repeat,  the  objective  is  to  generate  the  input  which  most  enhances  the 

1 

I  •  t  . 

estimates  of  and  Yp  from  the  input/output  data, 

I  ' 

|The  first  step  is  to  compute  the  elements  of  A.  and  B,  and  then  F. 

A  recursive  equation  for  w(i)  which  greatly  simplifies  the  calculation 

of  the  elements  of  Aj  and  B  is  derived  next.  By  the  Cayley-Hamilton  theorem, 

T  ’ 

4>  must  satisfy  its  own  qharacterjistic  equation,  so  it  follows  that 


-  ^jlp 


T,n-i2 


Postmultiply  by  h  and  premultiply  by  ($  )  yields 


,  T.n.  ,  T.n-1,  ,.T.n-2, 

(0  )  h,  -  -  02'^  ’  -  ■  ^  - 


(30) 


(31) 


It  follow's  from  the  definition  of  w(n)  that 

I  j 

w(n)  =  -  0pw(n  -  1)  -  0jw(n  -  2)  (32) 

for  n  2. 


The  second  step  is  to  compute  the  eigenvectors  of  F,  the  one 
corresponding  to  the  largest  eigenvalue  of  F  being  the  optimal  input. 


For  the  numerical  examples,  the  model  parameters  were  chosen 
to  simulate  an  underdamped  and  overdamped  second-order  system.  The 
optimal  input  sequences  were  computed  for  the  model  parameters  listed  in 
Table  I  for  N  =  20  using  a  standard  computer  subroutine  to  find  the  eigen¬ 
values  and  eigenvectors  of  F.  The  optimal  inputs  along  with  the  system's 
unit  response  are  graphed  against  the  index  i  in  Figs.  1  to  4. 


TABLE  I 

PARAMETERS  USED  IN  THE  EXAMPLES 


System 

Characterization 

h 

^2 

e 

Figure  No. 

-0.  50 

-0.  25 

0.  25 

0.  05 

1 

Overd.amped 

-0.  50 

-0.  25 

0.  50 

0.07 

2 

0.  50 

-0.  25 

0.  25 

0.  37 

3 

Underdamped 

0.  50 

-0.  25 

0.  50 

0.  21 

4 

The  following  discussion  describes  the  physical  significance  of  the 
optimal  input.  In  all  the  cases  where  ^2  ~  '  prediction  error  in  Xj(i)], 

the  majority  of  the  input  energy  is  spent  in  the  first  third  of  the  interval 
(Figs.  1  and  3).  This  is  apparently  an  effort  to  get  the  system  moving  as 
soon  as  possible.  In  cases  where  >  "  ^2  [prediction  error  in  Xj(i)  as 
well  as  X2(i)]  the  input  energy  is  more  dispersed  throughout  the  interval 
(Figs.  2  and  4).  This  dispersion  of  the  input  energy  appears  to  bo  a  hedge 


UNIT  RESP.  .  INPUT  US  INDEX 


UNIT  RESP 


Fig.  2.  Overdamped  with  k  = 


UNIT  RESP. .  input  us  index 


I'nclerdamoed  with 


Fig.  4.  Underdampcd  with  k  -  0.  50. 


against  the  poorer,  prediction  accuracy  of  the  plant  (reflected  in  the  larger 

I 

value  of  The  purely  random  input  is  better  with  faster  systems  than 

I 

with  slower  on.es,  but  the  effidiency  w^s  never  observed  to  be  over  .  50  in 
any  of  the  cases  tried.  ' 

The  majority  of  the  energy, of  the  input  appears  to  be  in  frequencies 
higher  than, that  contained  in  the  unit  response  of  the  system,  but  not  so  high 
as  to  be  greatly  attenuate<jl  by  the  system.  At  least  this  appears  to  be  one  of 
the  factors  taken  into  account  by  the  optimal  inpvit.  It  also  explains  why  the 
random  input  is  better  for  faster  systems  than  for  slower  ones:  more  of  the 

*  I 

energy  in  the  random  input  is  attenuated  bythe  .slower  system. 

The  analysis  and  the  numerical  results  described  in  this  example  arc 
useful  for  the  situation  where  one  has  fairly  accurate  prior  knowledge  of  the 
tr'ansienl  response  and  one -step-ahead  prpdictidn  accuracy  of  the  unforced 

.  '  I 

systems,  but  no  prior  knowledge  of  the  input  sensitivity  parameters.  In  this 
case,!  one  can  dptermine  an  input  test  signal  which  enhances  the  estimates  of 
the  input  sensitivity  parameters. 


I 


), 


I 
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VII.  INPUT  SYNTHESIS  FOR  THE  GENERAL  CASE 


In  general,  the  parameters  k,  and  ^  must  be  identified.  The 
residual  covariance  b  can  be  estimated  from  the  measurement  residuals 
after  the  identification  of  the  other  parameters  is  complete  [1].  Define^ 
as  the  3n  X  1  vector  of  unknown  model  parameters: 


e4 


The  elements  of  the  Fisher  information  matrix  in  the  general  case  are 


m. .  =  Tr 
1)  2 


92  _-i 
90.  90. 

1  J 


.  T  9A^  ^-1  9A 

-  90.  ^  90.  - 

1  J 


(33) 


The  expression  for  Tr[M]  follows  immediately: 


3n 

Tr[M]  =|]^Tr 
i=l 


,-l  92’ 


1  T 
+  f  uFu^ 
b  —  — 


(34) 


where 


9A 

90. 
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From  Eq.  (10)  it  follows  that 


which  shows  that  the  input  does  not  directly  affect  the  identification  of  k. 

Since  B  is  a  function  of  k,  the  optimal  input  will  depend  on  k,  but  it  is 
not  possible  to  enhance  the  identification  of  k  by  manipulating  the  input. 

This  makes  sense  because  it  is  well  known  that  the  Kalman  gain  is  independent 
of  the  input. 

The  problem  one  faces  in  maximizing  Tr[M]  with  respect  to  u  is  that 
8A/80.  is  a  function  of  ^  and  B  is  a  function  of  ^  and  k  which  are  unknown. 

One  way  to  synthesize  an  input  is  to  use  the  a  priori  values  of  ^  and  k  in 
the  optimization.  If  there  is  time  during  the  identification  process,  the 
input  can  be  reoptimized  using  the  partially  identified  values  of  ^  and  k  to 
replace  the  a  priori  values. 
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vm. 


CONCLUSION 


An  approach  to  optimal  input  synthesis  for  linear  system  identification 
has  been  presented.  The  significant  result  is:  the  optimal  input  is  an 
eigenvector  of  a  matrix  related  to  the  Fisher  information  matrix  and  in 
general  it  is  a  function  of  the  unidentified  system  parameter.,. 
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